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MODIFICATION OF A THREE-DIMENSIONAL SUPERSONIC NOZZLE ANALYSIS 
AND COMPARISON WITH EXPERIMENTAL DATA 


by Allan R. Bishop and Bobby W. Sanders 
Lewis Research Center 

SUMMARY 

A computer program previously developed to analyze three-dimensional supersonic 
nozzles by the method of characteristics has been modified to study less restrictive 
nozzle geometries and nonuniform inlet conditions. An example which is discussed in- 
dicates that averaging a nonuniform initial profile for one -dimensional calculation may 
generate significant errors. The analytical method and the data from a three- 
dimensional experiment are compared. The results show generally good agreement 
between analysis and experiment. 


INTRODUCTION 

There have been several methods developed to analyze three-dimensional super- 
sonic nozzle flows by the method of characteristics (refs. 1 to 4). Moretti (ref. 2), 
Rakich (ref. 3), and Dash and Del Guidice (ref. 4) use various approximations to the 
characteristic equations to produce a calculation accurate to first order in the grid spac- 
ing. The scheme of Ransom, Cline, Hoffman, and Thompson (ref. 1) is an exact solu- 
tion to the characteristic equations and is accurate to second order in the grid spacing. 
This report describes modifications to the computer code of reference 1 that allow it to 
accept a less restrictive wall geometry and initial profile. 

In a finite difference analysis it is often very convenient to specify the boundary con- 
ditions in an algebraic form. The values are usually continuous and well behaved, and 
differentiation and integration are straightforward. The computer program in refer- 
ence 1 allows the use of only a limited class of algebraic expressions for the solid- wall 
boundary; the axial variation in the throat is specified as a circular arc matched to a 
quadratic form for the expansion region of the nozzle. The cross sections are specified 
as superellipses, which may be different in each quadrant. 


There are, however, many nozzle geometries that are not easily specified in an 
algebraic form. Nozzles designed for optimum performance or corrected for boundary 
layer growth are usually described by a series of data points. Experimental nozzles 
are often designed for convenience in manufacturing rather than to fit any algebraic 
equations. Therefore, an alternate solid- boundary routine has been developed which 
retains the superelliptical cross sections but allows the axial geometry variation to be 
specified on a point-by-point basis. With this form of input a much more arbitrary 
geometry can be analyzed and complex nozzle geometries can be conveniently specified. 

At the initial station in the duct, an initial flow field must be specified. In refer- 
ence 1 an initial profile with radial variations in total pressure, total temperature, 

Mach number, and so forth, could be analyzed; but no azimuthal variation was permitted. 
This restriction has been relaxed to allow a general two-dimensional variation of flow 
properties at the inlet station. The initial flow field is specified on a point- by- point 
basis, and a two-dimensional interpolation scheme is used to determine intermediate 
points. While this procedure permits a very general initial flow field to be considered, 
some care must be taken to select conditions that are consistent with the boundary con- 
ditions. 

The first half of this report documents these modifications and gives the results for 
a few sample cases, while the second half describes a three-dimensional supersonic 
flow experiment and compares the experimental data with the three-dimensional analy- 
sis. 


PROGRAM MODIFICATIONS 
Changes in Geometry Specification 

The solid-boundary routines in the computer program of reference 1 have been 
modified to accept axial variations in the body coordinates specified on a point-by-point 
basis. The superelliptical cross sections have been retained, but the intercepts of the 
superellipses with the coordinate planes and the superelliptical exponents are specified 
as a series of points. The exponents need not be specified at the same locations as the 
intercepts. A variable centerline coordinate has been added, so that curved or S-shaped 
ducts may be studied. This variation is also specified as a series of points. 

Interpolation method . - Whenever a curve is specified by a series of points, an in- 
terpolation scheme is required to determine intermediate values. The method developed 
for the modified nozzle computer program uses Stirling's three- point interpolation 
twice: once with the center point behind the new point (as in region I of fig. 1) and again 
with the center point ahead of the new point (as in region II of fig. 1). A weighted aver- 
age is taken of the two results, the weighting factor being the distance along the abscissa 
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from the new point to the center point of the interpolation region. With this method of 
interpolation the function values and the first derivatives are continuous. This routine 
is particularly good near the junction of two regions with a discontinuity in curvature, 
such as a circular arc changing to a straight line. Splined interpolation will oscillate 
about the true value for some distance beyond the junction, while the method used herein 
provides a much smoother transition. 

Within the computer program the interpolation is performed in a separate subrou- 
tine and therefore may be easily replaced. Since the interpolation uses points on both 
sides of the unknown point, it is advisable to supply additional points on either side of 
the region of interest. While the subroutine will interpolate between the last two points 
and, in fact, will extrapolate beyond the last point, the accuracy is necessarily reduced 
in these regions. 

Example solutions . - To illustrate the new input scheme and to demonstrate its 
capabilities, two representative nozzles have been analyzed. A simple axisymmetric 
nozzle with a straight centerline is used for example I. This example serves as a check 
between the new and old input methods. A three-dimensional nozzle, demonstrating the 
flexibility of the new input, is used for example II. 

The nozzle used in example I is axisymmetric with a circular-arc throat and a 10° 
conical exit section. Figure 2 is a sketch of this nozzle, and the coordinates are given 
in table I. The uniform initial profile at the entrance station (throat) has a Mach number 
of 1.1. Figure 3 is a representation of the Mach number profile at the exit station. 
One-quarter of the exit plane is shown with the nozzle centerline at the origin of the 
coordinate system. The curves in figure 3 are elements of the surface formed when the 
Mach number is plotted as a function of position in the exit plane. This type of "sur- 
face” presentation is not necessary for axisymmetric flow, but it is very useful in vis- 
ualizing three-dimensional flow and is used in subsequent figures in this report. No 
differences between the results obtained from the old and new methods of specifying 
geometry were observed. 

Example II demonstrates the flexibility of the new input geometry. Figure 4 is a 
sketch of the nozzle used in example II, and the coordinates are given in table II. This 
nozzle has a varying cross section that changes from circular at the throat to almost 
rectangular at the exit. The variation is made by changing the superelliptic exponent of 
the z-coordinate from 2 at the throat to 20 at the exit. The exponent of the y-coordinate 
is constant at 2. In addition, the centerline of the nozzle has a slight S-shaped bend in 
the x-y plane. The uniform initial profile at the entrance station (throat) has a Mach 
number of 1. 02. Figure 5 shows the Mach number distribution at the exit station for 
example II. One-half of the exit plane is shown, and the nozzle centerline is offset from 
the origin of the coordinate system because of the bend of the nozzle. The asymmetry 
in the exit Mach number is caused primarily by the curvature in the nozzle centerline. 
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Initial Profile Modifications 


A three-dimensional calculation is necessary when an axisymmetric nozzle has a 
nonuniform inlet flow. The original code described in reference 1 would accept only a 
radial distortion of the inlet profile. This restriction has been relaxed to allow a gen- 
eral two-dimensional variation of flow properties at the inlet station. The initial pro- 
file is specified on a point- by- point basis, each point given by its coordinates and the 
flow properties at that point. There are no restrictions on the spacing or the distribu- 
tion of points, with the following two exceptions: all points must have the same axial 
coordinate and the points must be transformable into a rectangular grid. The points 
must be arranged in rows and columns corresponding to this grid, with every row having 
the same number of points and every column the same number of points. This latter 
restriction is a requirement of the interpolation scheme which is used to find intermedi- 
ate points on the initial- value surface. 

Two-dimensional interpolation method . - A two-dimensional interpolation scheme 
for the initial datum surface is necessary because the computer program chooses its 
own starting grid. The method used in the modified program determines the two trian- 
gular regions, I and II in figure 6, enclosing the new point (y^, Zj). A plane is fitted 
through the three corners of each region to determine the values of a flow variable as a 
function of y and z in each region. A typical flow variable S(y^, z^) is found on each 
plane, and the two values are averaged to find the final result. This method has been 
used successfully in several other problems. A continuous functional value is generated 
with a minimum of curve fitting and mathematical calculation. This interpolation is per- 
formed in a separate subroutine within the program, so that another interpolation rou- 
tine may be substituted. 

The flexibility used in the construction of the grid for the initial datum surface al- 
lows for a wide variety of initial inlet flow fields to be specified. However, care must 
be taken to match the flow field properties with the boundary conditions at the solid wall 
(velocity tangent to wall, etc. ). A provision is made within the program to alter the 
wall velocities slightly to match exactly the tangent condition at the solid boundary. 

Example solution . - To demonstrate the initial-profile scheme developed a third ex- 
ample is presented. Example III is an axisymmetric nozzle with a nonuniform inlet flow 
profile. Both Mach number and total pressure have a linear variation with the z- 
coordinate, giving the entrance flow both radial and circumferential distortion with sym- 
metry about both the y and z axes. Table HI gives the geometry of the nozzle, and 
table IV the nonuniform initial profile. A sketch of the nozzle used in example III is 
shown in figure 7. 

The Mach number distributions at the entrance and exit stations are shown in fig- 
ure 8. One- quarter of the exit plane is presented, with the nozzle centerline at the ori- 
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gin of the coordinate system. The distributions are not to the same geometric scale, 
since the exit radius is 4. 8 times as large as the entrance radius. The initial and final 
Mach number distributions are readily apparent in these two plots. 

A distorted Mach number profile at the exit station is not the only effect of initial- 
profile distortion. The performance parameters of the nozzle, such as mass flow and 
thrust, are also changed. In many situations involving airflow characteristics within 
nozzles, it is not convenient to consider local flow variations within the flow field. 
Therefore, the properties of the flow at one station are treated as though they are uni- 
formly distributed, so that one-dimensional equations can be used to predict the flow at 
other stations of interest. The use of averaged values, however, can produce signifi- 
cant errors in mass flow and momentum flux when compared to the actual integrated 
values. Table V shows such a comparison between the results of the three-dimensional 
analysis of example in and the results of a one-dimensional analysis that used several 
different averaging methods. A simple average, an area- weighted average, and a 
mass-flow-weighted average were used to determine a single value for Mach number 
and total pressure at the entrance station. There values were used in a one-dimensional 
analysis to find the mass flow and thrust. 

The three one- dimensional methods give different mass flows and thrusts and are in 
error from the three-dimensional values by varying amounts. * All three one- 
dimensional methods predict thrust values that are above the value for the three- 
dimensional calculation. This means that a one- dimensional calculation with an aver- 
aged profile at the initial station may give an optimistic value of thrust, and that any 
distortion of the initial profile may tend to degrade the performance of the nozzle. 


COMPARISON WITH EXPERIMENT 
Description of Experiment 

Some three-dimensional flow data have been obtained as part of an experiment in the 
10- by 10- Foot Supersonic Wind Tunnel at the Lewis Research Center. A series of 
plates were placed in the throat of the nozzle to generate distorted flow fields in the test 
section. A calibration rake was also installed in the test section, permitting a limited 
survey of the three-dimensional flow. Data from this experiment will be compared with 
the results of an analysis of the same geometry that were obtained by using the modified 
three-dimensional method-of-characteristics program. 

iThe simple average agrees best with the three-dimensional analysis for this case 
because of the initial profile considered. For an arbitrary initial profile a simple aver- 
age would be expected to be least accurate. 
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A sketch of the experimental geometry is shown in figure 9. The floor and ceiling 
of the nozzle were fixed, and the side walls were adjustable to produce various test- 
section Mach numbers. As indicated in figure 9, the disturbance generator was a flat 
plate placed on the floor of the tunnel at the throat. Plates of two different sizes were 
used with several different wall settings to vary the Mach number and the strength of the 
disturbance in the test section. Near the entrance to the test section a multiple -wedge 
calibration rake was used to survey the flow properties at several axial stations. 


Description of Analysis 

Several approximations were made in analyzing this geometry. The rectangular 
cross section of the actual nozzle was approximated by a superellipse with exponents of 
20. Since the program cannot analyze a separated region, a triangular-shaped solid 
boundary replaced the vertical flat plate (fig. 9). The height of the triangle was the 
same as that of the plate. A uniform initial profile with a Mach number of 1. 02 was 
assumed to exist at the plate location. The boundary layer and other dissipative effects 
were not considered in the analysis. 


RESULTS AND DISCUSSION 

Figure 10 presents the static pressure distribution along the tunnel ceiling when the 
sidewalls were positioned for a nominal test-section Mach number of 2. 4 with no plate 
present. The static pressure ratio (static to total) on the ceiling centerline is plotted as 
a function of the dimensionless axial distance from the throat (actual distance divided by 
tunnel height). Both the experimental data and the results of the analysis are plotted in 
figure 10 for comparison. 

The agreement between data and experiment is very good. The Mach number pro- 
files (not shown) are very uniform, the variation being less than 1 percent at any station. 
These results are also in good agreement with the experimental data. 

Figure 11 shows the disturbance generated when a plate 30. 5 centimeters high was 
placed on the floor at the throat. Again the static pressure ratio is plotted against axial 
location. The floor and ceiling centerline values predicted analytically, as well as the 
experimental data for the ceiling, are plotted. The strong compression created by the 
plate can be easily traced as it was reflected between the ceiling and the floor. The 
agreement between data and analysis is good except near the throat. The experimental 
data indicate that subsonic flow was present downstream of the geometric throat, which 
was at the axial location of the plate. Apparently, a vena contracta was present in this 
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region, so that the minimum flow area did not coincide with the minimum geometric 
area. 

Figure 12 is a plot of Mach number against axial location in the region surveyed by 
the calibration rake. The four plots represent the data from four wedges near the center 
of the tunnel. The vertical positions were at 30, 43, 50, and 70 percent of the tunnel 
height. The Mach number predicted by analysis is presented for comparison. Both the 
data and the analysis indicate a strong compression crossing this region. The compres- 
sion is shown more clearly in figure 13, which shows the predicted Mach number distri- 
bution at several axial locations in the region of figure 12. The comparison in figure 12 
indicates that the analysis predicted the magnitude of the compression fairly well but 
that the compression occurred over a larger axial distance in the analysis than in the ex- 
perimental data. This is probably the result of interpolation in the analysis, which tends 
to smooth any rapid changes. 

Figure 14 is a plot of static pressure ratio against axial location but with the side- 
walls positioned for a Mach number of 2. 0. The same 30. 5-centimeter-high plate used 
for the previous data (wall Mach number of 2.4) was placed at the throat. Because the 
sidewalls at the throat are further apart for a Mach number of 2. 0, a gap of 30. 1 centi- 
meters existed between each end of the plate and the sidewalls. Both the floor and ceil- 
ing centerline static pressure values from the analysis are plotted, along with the ex- 
perimental data on the ceiling. As with the Mach 2. 4 wall setting, the subsonic flow 
near the throat is apparently caused by a vena contracta downstream of the plate. The 
deviation of the experimental data from the analysis in the dimensionless axial region 
from 1 to 2 may result from local separation caused by the strong compression reflected 
from the ceiling. Kuehn (ref. 5) indicates that the predicted pressure rise occurring 
near an axial location of 2 is enough to cause the flow to separate. The experimental 
data also indicate that the last reflection of the compression wave from the ceiling oc- 
curred upstream of the predicted location. The reason for this is not clear; but the 
separated region previously mentioned, the gap between the plate and the sidewall, or 
the boundary layer effects neglected in the analysis are all possible explanations. 

The data in figure 15 are for a wall setting which gives a nominal test- section Mach 
number of 3. 1. A plate only half as high as previously used, 15.2 centimeters, was 
placed at the throat. The static pressure ratio is plotted as a function of dimensionless 
axial distance. In figure 15 the analysis is in better agreement with the data, particu- 
larly in the region near the throat. The cause may be the weaker compression which 
was present because of the smaller plate. 

The variations in Mach number at each axial station for the Mach 2. 0 and Mach 3. 1 
wall settings are similar to that presented for the Mach 2. 4 wall setting, except that the 
compressions occur at different axial locations, as indicated by the pressure distribu- 
tions. In general, the distributions predicted analytically agree with the results of the 
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experiment, except for a smoothing and stretching out of the compression in the analy- 
sis. 


SUMMARY OF RESULTS 

The modifications to the original three-dimensional method- of- characteristics pro- 
gram permit the analysis of a wider variety of nozzle geometries and nonuniform inlet 
conditions. The results for the sample case with a nonuniform inlet profile indicate 
that a one- dimensional calculation may have significant errors, no matter how the non- 
uniformity is averaged. 

The comparison between experimental data and the three-dimensional analysis 
shows good agreement, especially since the method of characteristics assumes an in- 
viscid flow field. While the analysis smooths out the sharp changes across a shock wave 
and ignores all dissipative effects, it does predict the location and magnitude of strong 
compressions fairly well. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, December 4, 1973, 

501-24. 
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TABLE I. - COORDINATES OF AXISYMMETRIC 


NOZZLE USED FOR EXAMPLE I 


Axial coordinate 
direction, 

X 

Coordinate perpendicular 
to axial coordinate, 

Y 

-0. 05234 

1. 00137 

-.03490 

1.00061 

-.01745 

1. 00015 

0 

1. 00000 

. 01745 

1. 00015 

. 03490 

1. 00061 

. 05234 

1.00137 

. 06976 

1. 00244 

. 08716 

1. 00381 

. 10453 

1.00548 

. 12187 

1. 00745 

. 13917 

1. 00973 

. 15643 

1. 01231 

. 17365 

1. 01519 

. 18000 

1. 01631 

.20000 

1. 01984 

. 30000 

1. 03747 

. 50000 

1. 07274 

1.00000 

1. 16090 

2. 00000 

1.33723 

4. 00000 

1.68988 
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TABLE II. - COORDINATES OF THREE-DIMENSIONAL NOZZLE USED FOR EXAMPLE II 
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TABLE III. - COORDINATES OF AXISYMMETRIC 


NOZZLE USED FOR EXAMPLE III 


Dimensionless axial dis- 
tance from throat, 

X 

Coordinate perpendicular 
to axial coordinate, 

Y 

-9.30000 

1. 81850 

-5.90000 

1.49250 

-2.95000 

1.31630 

. 00000 

1.25300 

2.95000 

1.31140 

6.00000 

1. 48700 

9.40000 

1.79790 

12.90000 

2.21750 

16.20000 

2.67970 

18.84999 

3. 08010 

21.39999 

3.47680 

23.95000 

3.87120 

26.50000 

4.24650 

29.25000 

4. 61110 

32.20000 

4. 94850 

35.50000 

5.26050 

38.89999 

5. 51030 

42.29999 

5.69070 

45. 70000 

5. 81270 

49. 09999 

5. 88920 

52.59999 

5.93330 

56. 09999 

5.95870 

59.99999 

5.97950 

63. 09999 

6. 00000 

64.00000 



66. 00000 



70. 00000 
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TABLE IV. - INITIAL- VALUE SURFACE DATA USED FOR EXAMPLE III 


| Coordinates perpendicular to axial coordinate Mach number Total pressure! 
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TABLE V. - THREE-DIMENSIONAL EFFECTS IN CALCULATING 


NOZZLE PERFORMANCE 


Type of calculation 

Entrance 

Mach 

number 

Entrance total 
pressure, 
N/m 2 

Mass flow, 
kg/sec 

Thrust, 

N 

Three-dimensional flow analysis 

Variable 

Variable 

4.72 

24.2 

Average of initial conditions, 
one- dimensional flow analysis 

1.31 

8. 32xl0 4 

4.72 

24.6 

Area- weighted average of initial j 
conditions, one- dimensional flow 
analysis 

1.35 

8. 61X10 4 

4. 81 

24. 8 

Mass-flow-weighted average of 
initial conditions, one-dimensional 
flow analysis 

1.35 

8. 63X10 4 

4. 80 

25.3 


y 



Figure 1. - Diagram of interpolation scheme used to find inter- 
mediate points on solid boundaries. 



Figure 2. - Sketch of axisymmetric nozzle used for example I. 
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Figure 3. - Mach number distribution at exit station for axisymmetric 
nozzle of example I. 


Exit station -7 



Nozzle 



Figure 7. - Sketch of axisymmetric nozzle used for example HI. 



(a) Entrance station. (b) Exit station. 

Figure 8. - Mach number distribution at entrance and exit stations for axisymmetric nozzle 
of example III. 
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Dimensionless axial distance from throat, X 

Figure 10. - Undisturbed static-pressure distribution in 10- by 10- 
Foot Supersonic Wind Tunnel. Wail setting, Mach 2.4. 



1 2 3 4 5 6 

Dimensionless axial distance from throat, X 


Figure 11. - Disturbed static-pressure distribution in 10- by 10-Foot 
Supersonic Wind Tunnel. Height of generator, 30. 5 centimeters; 
wall setting, Mach 2.4. 
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WcuQc 

location, 

percent 

— 70 
__ 50 

— 43 

— 30 


Tunnel cross section show- 



Figure 12. - Mach number as function of axial length for 
four wedges. Height of generator, 30. 5 centimeters; 
wall setting, Mach 2.4. 
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Ratio of local static pressure to tree-stream total pressure 



Dimensionless axial distance from throat, X 

Figure 14. - Disturbed static-pressure distribution in 10- by 10-Foot Supersonic 
Wind Tunnel for 30. 5-centimeter-high disturbance generator and wall setting 
of Mach 2. 0. 



Dimensionless axial distance from throat, X 

Figure 15. - Disturbed static-pressure distribution in 10- by 10-Foot 
Supersonic Wind Tunnel for 15. 2-centimeter-high disturbance 
generator and wall setting of Mach 3. 1. 
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